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SYNOPSIS 


In this study, hydraulic and thermosyphon models are 
formulated to investigate the transient and steady state flow 
and, thermal behaviour of a CAIPU Reactor core in the event 
of loss-of-flow accidents, resulting from pimiping failures. 

The hydraulic analysis leads to a flow coast down, which 
determines mass flow rate versus time. The transient t^pera- 
ture behaviour is evaluated nexb. In view of safety assessment 
following the transients, major attention is directed towards 
the performance of thermosyphoning (Natural convection) effects 
in the vicinity of a horizontal cylindrical fuel rod to extract 
the after heat due to the decay of fission products, under both 
isothermal and non-isothermal surface conditions on the fuel rod. 
It is found that natural convection Nusselt number is approxi- 
mately one-fifteenth of the forced convection (for coolant 
velocity 27.2 ft/sec) value. Furthermore, the surface tempera- 
ture variation has a significant effect on the heat transfer 
from the fuel rod. 



CHAPTER I 


INTRODUCTION 

1,1 SAFETY OF POWER REACTORS 

The thorough safety assessment of a nuclear power 
plant necessarily involves consideration of the potentiality 
and consequences of a large niomher of reactor accidents, 
ranging from minor mishaps to truely catastrophic events. In 
recent years the efforts associated with assuring the safety 
of large nuclear power reactors has accelerated rapidly. The 
specifics of the accidents sequences that must receive atten- 
tion depend strongly on reactors concept and on the details of 
design and operation of heat transport system, fission product 
barriers and safety systems. Among broad categories of reactor 
accidents, loss-of-flow accidents (LOPA) encompass one of the 
most conceivable nuclear power plant accidents, 

loss-of-flow accidents can result either from pumping 
failures or from obstruction in coolant loops. Flow obstruc- 
tions are most likely to lead to local flow starvation in the 
core, whereas pumping failures are most likely to affect a 
major part or all of the core. Pumps, however, are normally 
in active state of operation and they are subject to rapid 
failure either from mechanical seizure or from power failure. 

A catastrophic mechanical seizure of pumps may occur nearly 
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instantaneously, causing a precipitous drop of the pump head 
to zero. But, pump failure due to loss of power does not result 
in an instantaneous loss of head, since even if pump motor 
suddenly fails, the rotational inertia of the pump impeller and 
flywheel cause the pump to coastdown gradually over a period of 
time that is often of the order of several seconds. 

In loss-of-flow accidents, caused by a power failure 
to the all coolant pumps, the pressure drop across the core 
decreases, causing decreased mass flow rates and increased core 
temperature. Plow starvation may also take place, if there is 
an increase in the core hydraulic loss coefficient. Plow 
decrease results in deterioration of the heat transfer coeffi- 
cient and coolant flow rate. As a result, the elevation of 
the coolant temperature may lead to a boiling crisis and subse- 
quent insulation of the fuel element from the coolant. The 
equilibrium flow following the transient is quite small, since 
it is driven only by natural circulation. Thus, the most 
important action to be taken in the majority accident situations 
is the prompt shutdown of the nuclear chain reaction. This 
action rapidly reduces the power to the level of decay heat and 
therefore greatly alleviates the energy dissipation problem 
under adverse condition. 

Once reactor shutdown is achieved immediately, the 
primaiy problem is to assure an orderly path of removal of 
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decay heat. Otherwise, the energy imbalance would eventually 
lend to the meltdown of core and cause penetration of the 
primary system. Thus, a major engineered safety feature of any 
nuclear power plant is the emergency core cooling system. 

Equally important, in the event of a loss-of-flow accident is 
that, continuous core cooling must also be provided over a perioc 
of time required to actuate emergency core cooling system or in 
situations when emergency core cooling is not immediately 
available. The rate of decay heat production during this period 
after shutdown may still be capable of melting fuel elements. 

Under such circumstances, thermosyphon system is 
extremely important and capable of transmitting farely large 
rates of heat transfer in the earths gravitational field. It 
determines amount of core cooling .immediately available, even 
if the driving pressure following the transient fails to sus- 
tain natural circulation of coolant in the primary flow circuit. 
The reactor core consists of a number of flow channels, acts as 
a systan of thermosyphons and provides essential cooling 
required during such emergent situations. A study of the thermo- 
syphon systans reveals the nature of physical phenomena associa- 
ted with the process of heat transfer. 

1.2 THEBMOSYPHOIT SIBUm 

The heat transfer system which has become known as 
the thermosyphon has many applications in engineering, ranging 
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from its' use in Bakers' <ven to the cooling of turtine rotor 
blades, internal combustieri engines and nuclear reactors. A 
thermosyphon is a circulating fluid system, whose motions are 
caused by density differences in a body force field. Davis 
and Morris [1] have suggested that thermosyphon can be catego- 
rized according to (i) nature of boundaries (whether the system 
is open or closed to mass flow), (ii) the regime of heat transfer 
(whether the process is purely natural convection or mixed 
natural and forced convection) and (iii) the nature of the body 
force (gravitational or rotational). The thermosyphon system 
has an intrinsic function of removing heat from a prescribed 
source and transferring heat and mass over a specific path 
(frequently a recirculation flow) and rejecting heat and mass 
to a prescribed sink. These flows are intrinsically driven by 
thermal buoyancy forces either locally or in overall sense. A 
simple loop flow may be the result of local buoyancy forces, 
but a mult ibranched flow circuit may easily incorporate sections 
in which the flow direction is contrary to the local buoyancy 
forces resulting from the pressure created by the overall system 
buoyancy forces. 

The classification of thermosyphon systems mainly 

includess 

(1) A common single phase, natural convection open system in 
the form of a vertical tube open at the top and closed at the 
bottom. * 



(2) A simple single phase natural convection closed system in 
the form of a vertical tube closed at both ends. 

(5) fhe single phase closed loop thermosyphon. 

Vertical Open Thermosyphon i 

The vertical, open thermo syphon as shown in Fig. 1.1 
provides the starting point for considering thermosyphon systems, 
As can be noted that, the primary effect of heating the wall of 
an open thermosyphon should be to cause some type of flow upward 
along the wall due to the buoyancy effects and an associated 
return flow downward in the core via the law of continuity. It 
has been found that for large heat fluxes on the wall (as may be 
the case of nuclear reactors), the buoyancy forces are suffi- 
ciently intense near the wall so that a boundary layer regime 
is obtained. On the other hand for weaker heat fluxes, the 
buoyancy forces are less and the effect of shear is relatively 
enhanced owing to small inertial effects, causing the boundary . 
layer to try to fill the entire tube. Hence, compared to the 
ordinary boundary layer flow, the flow is impended. For suffi- 
ciently weak fluxes this effect is found to be rather uniform 
and a similarity flow is observed. The first experiments on an 
open vertical thermosyphon were reported by Foyle [2], but those 
were of limited extent. His observations show that the heat 
transfer was independent of tube length beyond a certain critical 
value although, it increases quickly with increasing diameter 
of the tube. He also pointed out that, for increasing values of 
length to diameter ratio (l/a), larger values of Rayleigh number 
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are required to attain any given heat flux level. 

Vertical Closed Thermosyphon ; 

The simple closed vertical thermosyphon as shown in 
Pig. 1.2 gained popularity as the problem of containment and 
pressurization became apparent in various application of open 
thermosyphon. It consists of a fluid completely enclosed in a 
tube situated in a body force field and lying parallel to the 
direction of the body force. Heat is applied to the ’lower' 
part of its length and cooled over the upper part, thus producing 
flox«r patterns as shown in Pig, 1.2. Refering to the figure, the 
first subscripts 0 and 1 associated with the temperature indica- 
te core region and wall respectively. The second subscripts 
1 and 2 refer to the top and bottom region of the thermosyphon. 

Q gives the heatflow. It has been found that, when closed 
themosyphon modelled carefully, it can be treated as a synthesis 
of two simple open thermosyphons approximately, joined at the 
mid tube exchange (coupling) region. Most of the modes of flow 
found for the open thermosyphon have been found for the closed 
thermosyphons. Primary problem concerned with the closed thermo- 
syphon is that of modelling the exchange region i.e. to find 
temperatures Tq^ and !Pq 2 » apply open thermosyphon results. 

Por lower values of temperature difference between the bottom 
and upper wall, the exchange mechanism was found to be basically 
convective and for higher values, the convective process is less 
stable and degenerates to a mixing exchange process. 
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Closed Loop Thermo syphon : 

A simple closed loop thermosyphon shown in Pig, 1.3 
is a similar geometric configuration, which can be found in many 
industrial situations. The virtue which made closed loop thermo- 
syphon interesting is that, it avoids chocking that occurs in 
open thermosyphon and diminishes exchange process as occurs in 
open and closed thermosyphons. The core boundary layer inter- 
action is also minimized. In principle, there is no limit to 
the types of flow that could be obtained in a closed loop 
thermosyphon, subject to various thermal, geometric and body 
force conditions. Under the condition of rotation of this 
configuration, for a given temperature difference between the 
wall and the coolant, gives a strong increase to the heat 
transfer coefficient. The better known application of closed 
loop thermosyphon, is cooling of internal combustion engines 
and nuclear reactors, when driving pressure establishes a 
natural circulation. 

1*5 BMim OP LITERATURB 

The veirtical thermosyphon systems and its various 
applications in engineering has been the subject of a good deal 
of research in recent years, although, probably because of the 
complexities of thermal buoyancy forces and their dependence on 
geometrical configurations and flow regimes, little attention 
has been given to it in its horizontal form. More attention has 
been directed towards the open thermosyphon in which the heat^ed 
lower region opens at the top into a large reservoir of fluid 
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from -which, heat is removed by allowing evaporation to occur. , 
Cohen and Baley [3]» reported a number of such experiments 
using both static and rotating apparatusy concentrating mainly 
upon the evaporative system. The original application of this 
system to turbine cooling was given by Schmidt [4]. Experimental 
investigation to remove heat from the reservoir by means of an 
immersed cooling coil was carried out by Cohen and Hartnett 
et al. [5>6]. Under these conditions, it was found that, the 
iluid returning from the reservoir has free access to the "core** 
of the thermosyphon. 

The analytical study presented by Light hill [7] has 
served as the foundation for most vertical thermosyphon analyses 
(laminar) upto the present time. His treatment of flow in the 
open thermosyphon yielded three laminar flow regimes, depending 
upon wall heat fluxes. He considered basic laminar boundary 
layer flow equations for an open circular thermosyphon with 
variable fluid properties as follows s 
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where, u, v are velocities in x and R directions respectively, 

g is gravitational acceleration. K, P, p, C M’ are 

thermal conductivity, pressure, coefficient of thermal expansion, 

density, specific heat and dynamic viscosity respectively. 

Subscript 1 refers to the bottom of the thermosyphon. 

With these equations, he treated the case of Pr = ...Vi, 
and thus neglecting the inertial terms after non-dim ensionali zing 
the equations and estimated that, the analyses would give only 
10 y, error for Pr = 2. For different values of wall heat fluxes 
at the bottom of the thermosyphon, he used Blasius-Howarth 
Series solution to solve the problem. He finally established 
theoritically the existence of boundary layer flow regimes for 
large heat fluxes, impe:^ded flow regimes for weaker heat fluxes 
and similarity flow regimes for sufficiently weak fltixes, 

Martins' [8] experimental results provided an excellent verifi- 
cation of lighthills' flow models. His local heat transfer 
results for lower values of wall heat flux also verified the 
existence of laminar impended regimes and similarity regimes. 

Light hill [7] again considered that the closed vertical 
thermosyphon system could be treated as a synthesis of open 
systems coupled together hydrodynamic ally in some manner by a 
region of limited extent. He suggested that complete mixing of 
the opposite streams in the coupling region would probably be a 
reasonable approximation. But it was left to Bayley and Lock [9]> 
who made the first comprehensive experimental and analytical 
study of the closed thermosyphon. I’heir attempt was to analyse 
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the system for the laminar houndary layer regime. They found 
that, in each section of the tube a thin laminar boundary layer 
was formed adjacent to the tube wall, thus producing a central 
core of fluid which moves slowly in the opposite direction. 

Three limiting possibilities were suggested for the manner in 
which the open systems were coupled together depending upon the 
fluid properties. Those were ; 

(1) The violent collision of the axially converging boundary 
layers giving rise to vigorous mixing, 

(2) The non-int eract ive crossing over of the boundary layers, 
each becoming the core of the opposite section. 

(5) The return of each boundary layer to its own section; 
thermal energy is thus transferred across the coupling 
region by axial conduction only. 

They theoretically treated the above mentioned hydro- 
dynamic couplings of the mid tube exchange mechanism and deter--, 
mined the effect of Prandt|Ll number on heat transfer rates. 
Considering the integral form of the energy and momentum equation 
as cited below, 

1 1 

^00 = f (T-«i)d{y/6) - fe ® (1.3.6) 

1 

/ ’i(T-Ii)4(y/6) (1.3.7) 

O 

where, 6 is the boundary layer thickness. Other notations have 
their usual meaning and subscripts ’o’ and '1’ refer to wall and 
core respectively. 



The approach developed hy them permitted solutions to 
he obtained for vertical tubes of almost any cross section and 
including all fluid for which 1 <'r Pr The method is based 

on the assumptions that, the two dimensional form of the boundary 
layer equations are approximately vc-lid near the tube wall and 
using Vonk^^arman - Pohlhausen technique, they obtained the 
following expressions. 


- X / I 1 

( = 0.677 (pr + 20/21^ ^ Pr(t 


for an isothermal wall and. 
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Pr(t,)„ 1/4 I 

"" ^’^^^Cpr + 4/5^ ^ ~ 
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[■ 


0 . 410 ( 1 ^ - /^^ ^) - 2 j 


(1.3.9) 


for a uniformay heated wall. Subscripts 'd' and 'H' refer to 
the diameter and heated section respectively and t^ is the modi- 
fied Rayleigh number. 


Japikse and Winter [10] pointed out that, the wall 
temperature on each tube half should be used to model the flow 
processes in that tube half. This is of considerable importance 
to avoid error in heat transfer rates, in situations where non- 
isothermal wall condition exists. 

Howevefir, the application of the above studies, mostly 
dealing with "vertical" thermosyphons are basically ccnfined only 
to light Water Reactors (LWR). Boiling Water Reactors (BWR) a 
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specific case of L¥Rs, have a vertical core conf iguration and 
easily actuate vertical thermosyphon cooling processes in the 
event of flow failure accidents. On the other hand, thermosy- 
phoning for fluid enclosed in a horizontal tube or channels 
situated in a body force field and lying perpendicular to the 
direction of the body force on a horizontal plane, is known very 
little. Thus, for "horizontal" reactor core configuration e.g. 
CAiroU Reactor (RAPP), this aspect is of considerable importance 
to assess the safety features involved in thermosyphon cooling 
in flow failure accidents. 

No literature, however, found available, explicitly 
treating horizontal thermosyphon systems in the gravitational 
force field. Since, thermosyphoning is closely associated with 
natural convective heat transfer processes, some insight can be 
gained by reviewing natural convection processes around the 
horizontal objects. 

Por laminar flow around horizontal cyrlinders, Hyman 
et al. [ll] found that the following equation for the average 
Nusselt number correlated data well for water, tolune, silicone, 
mercury etc. 

NU[j = 0.53[( p^ /o.952 ^ (Srp-Pr)]^/'^ (1-5.10) 

In some cases of horizontal cylinder, the entire surface may not 
be parallel to the body force caused by gravity. The variations 
of the angle between the surface and gravity vector was considered 
by Herman [12], who derived an equation for the local Nusselt 
number at various positions around a horizontal cylinder for air 



(Pr = 0.74). The equation is 

KUp = 0.604 0(9) (1.3.11) 

The function 0(9) depends upon the value of 9, measured from 
the bottom of the cylinder. An integral approach to this problem 
by levy [13] and a similarity solution approach by Eckert and 
Soehngen [14] > also agree well with the calculated values by 
Herman. 

On the basis of the foregoing reviews, the prime 
objective in the present investigation is to determine, thermo- 
syphoning effects, to provide immediate core cooling in Reactors 
of horizontal core systems in the event of loss-of-flow accidents. 

A typical example is the CAHliU (HAP?) type Reactor. 

1.4 DESCRIPTION OB' CAEDU REACTOR gfSTEM 

CANDU is a heavy water moderated pressurized heavy 
water cooled, U 02 fuelled pressure tube reactor. The reactor 
has a horizontal cylindrical tank and tube assembly. This assem- 
bly is known as calandria. Natural Uranium in the form of 
pressed pellets is sheathed in Zircaloy tubes to form cylindrical 
fuel elements. Nineteen such elements are assembled to form a 
fuel bmdle. Calandria tube contains thick walled coolant tube 
having a fuel bundle and primary coolant i.e., heavy water at 
high pressure. The calandria shell contains heavy water at hi^ 
pressure and it is used aS a moderator. 

Heavy water essentially serves two purposes, namely 
as a moderator and as a coolant, each having its own closed 
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circulating system. The fuel coolant system is called primary 
heat transport system is a high pressure and high temperature 
circuit with no bulk boiling. But, the moderator circuit is a 
low temperature low pressure circuit. Both circuits are insula- 
ted from each other in the reactor core by a gas annulas between 
the calandria tube and coolant tube. 

Heat generated in the fuel elements is transferred to 
heat exchangers in the primary heat transport system, where the 
heat is utilized to generate steam in a secondary coolant system, 
consists essentially of two identical loops in series, one at 
each end of the reactor, I'our boiler assemblies and four pumps 
in parallel are included in each loop. The primary coolant is 
circulated through alternate coolant tube assemblies in opposite 
directions. Each boiler ass^bly consists of ten vertical 
U type tubes and shell heat exchangers in parallel and a steam 
drum. The hot heavy water from the reactor enters the tube 
sides of the heat exchangers at the bottom of their recircula- 
ting sections. Then it returns from the bottom of the preheater 
legs of the heat exchangers to the reactor. Schematic diagram 
of the primary flow circuit is shown in Pig. 1 , 4 . 

The secondary coolant system consists of boiler feed- 
water entering the shell side of each heat exchanger at the hot ton 
of the preheat leg. The feedwater is heated near to the satura- 
tion temperature. As a result, it boils in the upper portion of 
the leg and, enters the steam drum above it. Water and steam are 
separated here. The water at saturation temperature flows through 
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downcomers to the bottom of heat exchanger recirculating leg, 
where it boils and returns to the steam drum. Saturated steam 
is supplied to the turbine unit, hiecharge valves are provided 
in the steam mains to release steam to dump condenser in the 
case of pressure rise in steam drum. Such a facility is provi™ 
dedj in the case turbine is unable to accept enough steam. The 
turbine and reactor have adequate controls for safety and smooth 
operation of the power plant. 

1 • 5 SCOPE OF THE THESIS 

Scope of the present work is mainly centered on the 
investigation of thermo syphoning effects in CAIEU Reactor System. 
It determines amount of core cooling immediately available as an 
essential safety requirement in loss-of-flow accidents particu™ 
larly, in situations, when emergency core cooling is not imme- 
diately available and equilibrium driving pressure ceases to 
establish a natural circulation in the primary circuit. However, 
transient flow and temperature behaviours are also essential- 
features of a loss-of-flow accident and must be evaluated prior 
to steady state analysis. 

Under operating conditions the bulk of the fission 
product inventory remains trapped in the fuel. The fission 
products that escape from the fuel consists primarily of fraction 
of noble gases and other species of more volatile fission product 
These collect in the fuel-cladding gap and in other void regions 
on the cladding envelope. Their existence may affect the gap 
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conductance and give rise to non-isothermal cladding surface 
conditions, resulting in a variation of steady state heat 
transfer across the surface. 

In the investigation proposed here, the objectives 
will thus be to ; 

(a) Study the flow coastdown during the transient, which 
determines decrease in mass flow rate with time. 

(b) Pind out whether the driving pressure, following the total 
loss of pumping power can establish a natural circulation 
equilibrium flow in the primary circuit. 

(c) Determine the temperature distribution of the core during 
the transient and approximate evaluation of the time 
required to attain steady state. 

(d) Evaluate the important thennosyphoning (steady state natural 
convection) effect on the flow and heat transfer in terms 

of the distributions of velocity and temperature in the 
adjacent boundary layer regime of fuel elements, considering 
both isothermal and non-isothermal cladding surface conditi- 
ons. Cases of linear and non-linear variations of surface 
temperatures are treated in non-isothennal conditions. 

Obviously, detailed analysis of a project of this 
magnitude would be a lengthy one and requiring much time and 
effort, I'he analysis will thus be carried out only as far as 
is required to arrive at a qualitative conclusion for the objec- 
tives listed in items from (a) to (c). In the perspective of 
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thermosyphoning, listed in item (d) , an asymptotic similarity 
solution of the equations of continuity, conservation of momentun 
and energy in the boundary layer regime will be sought by 
numerical methods f or an axisymmetric, constant property, laminar 
boundary layer flow of a viscous, incompressible and newtonian 
fluid aro\md a horizontal cylindrical fuel element. There is an 
arbitrary distribution of temperature on the surface of the fuel 
rod. Gravitational force field generates buoyancy forces 
responsible for thermosyphoning in the adjacent boundary layer. 

For ease in formulation and analysis of the present 
problem as envisaged in the proceeding objectives, those may be 
mainly divided into two parts. The first part named as ” Flow 
modelling’, comprises first two objectives (a) and (b), which 
deals with transient and steady state flow behaviours in the 
primaiy system and the second part namely, Thermal modelling*' 
comprises rest two objectives (c) and (d) , deals with transient 
temperature distribution in the core followed by thermosyphon 
cooling. 
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CHAPTER II 

PLOW MODELLIHG 


The determination of both transient and steady state 
flow behaviour following the transient, associated with reactor 
accidents, due to common mode power failure to pumps requires 
formulation of a hydraulic model for the reactor core and primary 
circuit. The mass flow rate versus time must be determined to 
gain an insight into the flow coastdown phenomena. Eollowing 
the transient, equilibrixun flow is quite small and driven only 
by natural circulation. It is required to verify that whether 
natural circulation can be sustained in the primary circuit. 


2.1 HYDRAULIC TRANSIENT MALI SIS 

To accomplish objective (a) of Sec. 1,5> transient 
hydraulic analysis is formulated. Refering to Bird et al. [15] 
and lewis [16] j control volume form of the continuity and 
mechanical energy balance equations in Eulerian frame are utili- 
zed to arrive at a relationship between mass flow rate and 
pressure drop for a flow channel. 
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[fj + [f 

c^v* c«s« 

• VP dV + jjj v-g'dV - E^+Ep (2.1.2,) 

C.V. G.V, 

where, ^ = fluid density, V = fluid volume, s = surface of the 

A — > 

volume V, n = outward normal to s, v = fluid speed, v = fluid 

velocity, P = pressure, g = gravitational force vector, 

E^ = Rate of mechanical energy dissipation due to viscous force 

and E = Rate of mechanical energy generation due to pump. 



Pig. 2,1 ; Plow channel. 

Let us consider a specific flow channel as shown in 
Pig, (2.1),e is the distance along the centreline and dV=:A(e)de, 
where A(e) is channel area at e. Pluid velocity, speed, density 
and pressirce are considered function of e and time. 

V = v(e,t), V = v(e,t)e 

/?= f(e,t), P = P(e,t) 

Considering the volume V between an inlet and outlet e^, 
equations (2,1.1) and (2.1,2) may be written as 




.d 


fO 

j ^ /^'(e,t)v^(e,t)A{e)de + J f(eQ,t)v(e^,t)\(e) 

= - J v(e,t)A{e) de - g I /^(e,t)v^(e,t)A(e)de~E[^+E^ 






(2.1.4) 


where, g is the magnitude of gravity vector and v vertical 

2j 


component of v. 


In the absence of boiling, we can treat primary system 
flow transients by assimiing incompressible flow, i.e. 
which means that f> = f{z). Time independence of f causes the 
^irst term of equation (2.1.3) to vanish, and we have 


= /’(ej^)v(e^,t)A(s^) (2.1.5) 

Siiice, e_ is arbitrary, it gives mass flow rate /2(e)v(e,t)A( e)tc 
independexit of e. Thus, mass flox'f rate at any point along 
the channel is given as 


W(t) =; f{z) v(e,t) A(e) 


( 2 , 1 . 6 ) 


Using these incompressible flow approximations and 
carrying out the mathematics involved, the expression for pressm 
drop follows from equation (2,1,4) 




where, p « I’i”I’o» ^ = ^i“®o» P 


■ 1 

-u. 


(2,1.7) 


f (e)A(£)de, 
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1 _ r f dg 

= - L f^TA(e) » 



H 



9 


is an empirical loss coefficient and subscripts ’o’ and 'i' 
refer to outlet and inlet respectively. 


The sum of contributions to the pressure drop come 
from inertial, 'Viscous, hydrostatic and pump energy/ term 
respectively. Combining, viscous and acceleration pressure 
drop expressions, eqn. (2.1,7) can be written in a contracted 
f orm. 




k 

I 




2fk 

11-)^- (li- ^ 


( 2 . 1 . 8 ) 


where, K = e-, + ( « . 

r o'^o 

of the channel, and Az = 


Pi'^i 


) , known as hj/draulic resistance 


z “ z . . 

o i 


2.2 TREATMENT OF THE PRIMASY PLOW CIRCUIT 

The standard pressure drop and flow relationship given 
in equation (2,1.8) can be applied to CANDU system. Considering 
a simplified assumption that all pumps behave identically 
following a power failure and pump inertia following the loss of 
power is taken into account in the expression for pressure drop. 

With the above assumption, the accident can be treated 
in terms of the reactor core and N identical coolant loops (The 
loop consists of series of paths through the hot and cold legs and 
the heat exchangers). Core is characterized by a loss coefficient 
and flow area A^ . Each of the coolant loops has a loss 
coefficient and flow area A^i Normalized flow resistance i.e. , 
^c ~ ^1 ~ ^ defined to eliminate flow area 
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from the pressure drop equation. 

Excluding hydrostatic term for a horizontal core 
^z = = 0) in the pressure drop relationship, we 

write core pressure drop as 




j. fr 

^A^c dt ^c 2^ 


( 2 . 2 . 1 ) 


and for coolant loops 


Z^P-, = 


(t) 

dt 


E] -4r + g(fAz), 
2p ^ 



( 2 . 2 . 2 ) 


Where, W^ and W^ si’s flow rates in the core and in each flow 
loop respectively. The term (”—) multiplies the initial pump 

Qq 

head in order to approximate the decreasing speed of the pump 
during the transient. is the initial pimp speed and O is the 

decreasing speed during the transient. 

Applying, Kirch offs laws for pressure drop around any 
loop of the channel and hydraulic junction law for flows in a 
hydraulic circuit , we have s 


AP^ + = 0 (2.2.3) 

¥ = N¥^ (2.2.4) 

Now, combining expressions (2.2.1) and (2.2»2), we find that 
the core flow is determined from 

where for the primary system. 



(k) _ (Si) ,1 /L) 


K = K 


+ IC, and if' ^z) = {f /2>z)-, 

•RT ^ j- U -t. 


U 


SuTd script p refers to primary circuit. 


Hydrostatic head is small compared to pump head during 
the first stage of the transient and can be neglected from 
equation (2,2.5). It is only predominant after pump head vani- 
shes at a later stage. 


Refering to Burgreen [17] to model the pump speed, 
we write the equation of motion s 


I = - cCj2(t) (2.2,6) 

subject to initial condition <^(t) -(Jo ^ 

moment of inertia of the pxamp, c is a loss coefficient relating 

to windage torque. Solution of (2.3.6) gives 


CJ = 


Cdo 


1 + t/t. 


(2.2.7) 


Where, t^ = l/c Uj^ is the pimp halftime. Substituting equation 
(2.2,7) in (5,2.5) and neglecting hydrostatic term in accordance 
with previous argument, we obtain 


(t) 


d¥ 


+ K 


^ 11 = 


Pg H, 




A^p dt P 2 P (1 + t/t )2 


( 2 . 2 . 8 ) 


subject to initial condition 


¥ 


K H 


0 


(2.2.9) 
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The nonlinear differential equation can he written in terms 
of the normalized core flow utilizing the initial condition, as 


t (^J-) , _ 1 

1 dt ^ )2 


( 2 . 2 . 10 ) 


wh er e , t , = 


2 F 


WoKp(p. 


P 


designated as loop halftime. 


A solution of equation (2.2.10) gives the flow coast- 
do-^m relation in terms of the normalized flow as a f\mction of 
time, for CAMDU reactors designed with sufficiently large 
flywheels, the pump inertia may he made much largers than that 
of the fluid, and hence this is the case, the 

normalized flow rate varies much slowly with time and it is 
reasonable to neglect the time derivative in equation (2.2.10). 
This approximately gives flow coastdown hy 


W. 

W, 


1 

1 + "t7t 


( 2 . 2 . 11 ) 


For a typical recirculation pump, the pimap halftime, t is 
calculated to he 5 seconds. Table I gives numerical results 
of variation of normalized flow with respect to time, normalized 
with pump halftime, using relation (2.2.11). A flow coastdown 
plot of these results Is given in Figure (2.2). 

I'/hen t^ and t^ are of the same order of magnitude, 
a solution of nonlinear equation (2,2.10) must he sought. As 
such, no standard method in terms of elementary functions found 
applicable to it. Some aspects of the solution is given in 
Appendix A. 







Table I 


( Plow Coastdown Resvilts ) 


Pump Halftime 



5 sec. 


Ratio of the elapsed 
time to the pump 
halftime, (t/tp) 

j Pr auction of initial 
1 flow, ¥/¥o 

1 (Normalized flow) 

0 

1.0 

1 

0.5 

2 

0,33 

3. 

0.25 

4 

0.20 ' 

5 

0.16 

6 

0.14 

7 

0.12 

8 

0.11 

9 

0.10 

10 

0.09 

11 

1 0.08 

12 

1 0.076 

13 

0.071 

14 

: 0,066 

15 

0.062 


2.3 DRiyiH& PRESSmiE IH HATUSAL CIRGULATIOH 

The second important aspect in flow modelling i.e, , 
objective (b) Sec. 1.5» is to verify whether natural circula-tion 
can be sustained in the primary circuit following the total loss 
of pumping power. Because of the density differential between 
the coolant in the downcomer and that in the riser, a driving 
pressure is established. Provided the reactor is shutdown fast 
enough’ to prevent any damage, a central question thm becomes 
whether an equilibrium state can be established for natural 
circulation to remove decay heat. 
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To establish- such an equilibrium state the driving 
pressure must equal to the total system losses at the desired 
rate of flow of coolant in the channel. This must be verified. 
Following El-Wakil [18], the system losses are given by 




(2.3.1) 


Where, = Total pressure losses in Ibf/ft^, 

^ AP^= Sum of frictional pressure losses, in core 
riser and downcomer, 

y = Sm of the acceleration pressure losses, 

= Sum of the pressure losses due to resistance 
to flow such as, due to drag of submerged bodies- 
like spacers, fuel element handle etc. 

The total system losses found out from RAPP data is 
1.413 Ibf/in . The driving pressure may be evaluated as followss 


Z^P^ = (hydrostatic pressure in downcomer) 
- (hydrostatic pressure in riser) 



( 2 . 3 . 2 ) 


where, H is the hei^t of the downcomer as well as the riser 
and^ are mean densities in the downcomer and riser res- 

pectively. g is the gravitational acceleration and g_ is the 
conversion constant. 
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In standard gravity g/g^ = 1.0, Substituting 

V-* 

corresponding values in, equation (2.3.2), AP^ is found out 
to be 26.31 Ibf/ft^ i.e. , 0.182 Ibf/in^. 

Now, in the case the computed driving pressure being 
less than the total system losses, the coolant rate of flow 
decreases and natural circulation cannot be sustained. The 
flow xiLtimately comes to a stagnation as the steady state is 
attained. 
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CHAPIBR III 
THEMAL MODELLING 

A comprehensive analysis of the problem of thermal 
modelling for a reactor core under accident conditions and 
steady state immediately following the transient, frequently 
requires the numerical solution of large systems of partial 
differential equations for transient heat conduction, fluid flow, 
transient and steady state convection and a host of other 
phenomena likely to occur during and after the course of a 
transient. Simple lumped parameter models basing on energy 
balances, however, often provides sufficient accuracy when the 
primary objectives are to gain an understanding of the inter- 
action of various phenomena involved and to make rough estimates 
of the consequences of the accident. But, the present formula- 
tion of the problem largely revolves around two important aspects, 
namely, the transient temperature distributions in the core and 
steady state thermo syphoning cooling effects following the 
transient, which is of paramount importance inview of safety 
assessment . 

3.1 TRANSIENT gEMPERAlURE DISTRIBUTIOM 

Accurate description of the temperature distribution 
that appear in the reactor core under accident conditions 
requires that the time dependent mass, momentum and energy 
balance conservation equations to such situations, Mayer [20], 
must be solved successively by . numerical methods. The complexity 
in the process of solution is due to the time involvement and 
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severe nonlinearity for the presence of convective 
terms in the inertial force expression, This can be further 
agipravated owing to the reactivity effects throu^ tonpera- 
ture feedback. 

For the present investigation as proposed in the 
objective (c) sec. 1.5, we try; to formulate a simple trans- 
ient conduction analysis to arrive at an expression for the 
solution of this equation subject to given boundary 
conditionsi In the process, we shall deal with the analytical 
method to derive a solution, in order to gain an insist 
into the procedure of . solution for non-horn ogeneous boimdary 
conditions, limited resources of data and inadequate 
informations regarding the boundary conditions of the reactoi 
core during the process of transient, puts a constraint for 
arriving at a ccmplete solution of the problem. The analy- 
tical treatment will be general, and a solution can be 
immediately sought, once the required boundary conditions 
are known. 

In the case of veiy severe transients, the 
analytical technique may not provide the desired accuracy 
due to large changes in fuel tomperature and coolant 
conditions. This fact coupled with obvious compleEities of 
analytical technique and availability of digital computers 
has led to an increased use of numerical techniques. The 
math^aatical apparatus required for the finite difference 
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soliition of parabolic equation of heat conduction type 
and required convergence criteria for explicit finite 
difference solution was comprehensively treated by Crank 
and Nicholson [21], 

In treating the problem of transient temperature 
distribution in the- reactor core by a conduction analysis, 
some simplifying assumptions are made. 

These are s 

1. The radial dimension of the fuel bundle placed in the 
core is small compared to the axial i. e. horizontal 
dimension. 

2. V/hile the axial heat flux distribution (i.e. volumetric 
thermal energy) in the core varies sinusoidally, the 
radial distribution shows relatively small changes. 

Thus, the neutron flux and consequently the volumetric 
thermal source strength q*” , are considered as functions 
of only the axial coordinate. 

3. It will be also assumed that the physical properties 
of the coolant and the fuel (e.g. density, specific 
heat, thermal conductivity) remain constant and 
independent of temperature. 

4. Core is considered to be of constant cross-sectional 
dimension. 


I 
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The sinusoidal core heat flux distribution can 
be utilized to find out the initial temperature distribu- 
tion in the core with the help of the bomdary conditions 
at the inlet and outlet prior to the transient. The volu- 
metric thermal energy distribution is shown in Fig. 3»1» 

The distribution is of sine form with respect to the origin 
situated at the inlet and positive axial co-ordinate is 
taken towards the outlet. T^^ siid T 2 sere the inlet and 
outlet temperatures respectively, prior to the transient. 
Their numerical values are quoted in Appendix III. These 
values refer to the reactor conditions during steady state 
operation i.e. prior to the initiation of transients. 


x=o 



Inlet_^ ^1 
Fig, 3* Is 




(x) 






Extrapolated 

Curve 


L—Ijs. •X= & 

CORE To Outlet 

£ 

Sinusoidal Heat flux distribution in the core. 


The total length of the core is considered to 
be i? - 16.404 ft. (This includes extrapolation distance.) 
Under transient conditions with an initial sinusoidal 
volumetric thermal energy distribution in the axial direc- 
tion, the general conduction equation for heat flow 
considered only in one dimension (i.e. one space' variable) 
assumes the well Isnown form. 
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(3.1.1) 

p at 

where, t represents time, x axial coordinate, C specific 

P 

heat, /^density and Z thermal conductivity and T(x,t) is 
the tGmperatm*G distribution. The initial volumetric 
thermal energy distribution assumes the following form 

fit 

M = <l" sin (^) (3.1.2) 

where, q”' is the peak value of the thermal energy inside 
the core. 

The initial and boundary conditions to this 
boundary value problem is stated as followss 

T (0 , 0) == 

at t = 0 (3.1.3) 

T ( ^, 0) = T2 
T (0, t) = aj_(t) 

for t > 0 (3.1.4) 

H iJ0, t) - a2(t) 

A complete solution of the present formulation to 
evaluate the transient temperature distribution i.e. seeking 
a particular solution to the conduction equation demands 
that the- £mctions a^(t) and a 2 (t) should be explicitly 
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defined as a function of time. Ho experimental data or 
theoretical estimates of the natijre of functions a^Ct) 
and ttgCt) during the process of transient is available. 
However, a general solution to the problem can be sought 
by analytical method such that the transient tonperature 
behaviour during the course of flow transient can be known, 
if the functions specified. 

Initial temperatu re d istribution in the c ore ; 

The initial temperature distribution can be 
evaluated from the steady state conduction equation and 
thermal energy distribution, which represents a source 
factor. 

■ Thus, from equation (5.1.1) and (5.1.2) for 
t = 0, we have 

sin (^) = 0 (3.1.5) 

dx ^ ^ 

Above equation integrates twice to give 

T(x,o) = sin (^) (3-1.6) 

where, and C 2 are constants of integration. 

Initial taiperature distribution at any position 
can be determined by evaluating these integration constants 
with the help of the conditions enumerated in eqn. (3.1.3). 
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This gives two simiiltaneous equations in two unknowns 
and C 



which is solved to give, 

C 3 _ = T 3 _ and 

Hence, from equation (3,1.6), the core initial temperature 
distribution follows as, 

T(x,o) = )x + (|)^ sin(^) (3.1.7) 

Analytical method ; 

The complete boundary value problem will be: 

_ _K_ d^T(K.1:) 

'3* f°p 9x2 

With boundary conditions; and initial conditions ; 

T(x,o) = T^ + + (|)^ sinC^) = 

T(o,t) = a^(t), T(^, t) = a^(t) 

For non-homogeneous boundary conditions prescribed in the 
given problem, it becomes necessary to split up the boundary 
conditions into groups of partial conditions, and solutions 
of the differential equations satisfying each group separatel 
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must be sought. Then, their sum automatically satisfies 
the heat conduction equation, on accounts of it*s linearity. 
The original division into partial conditions should be 
carried out in a correct manner such that, the boundary cond 
tions will also be satisfied. 

The scheme for the decomposition of the boundary 
conditions and of the complete solution T(x,t) in partial 
solutions v(x,t) and w(x,t) is carried out as follows. 

T(x,t) = v(x,t) + w(x,t) (3.1.8) 


Then equation (3.1.1) becomes 

The transformed boundary conditions will be : 
v(o,t) + w(o,t) = a 2 (t) 
v(£,t) + w(^, t) = a 2 (t) 


(5.1.9) 


(3.1.10) 


and the initial condition becomes, 
v(x,o) + w(x,o) = T(x) 


(3.1.11) 


Now, w(x,t) can be chosen such that, 



= 0 


l^x^ 


w(o,t) 

= a^(t) 

(3.1.12) 

w(^,t) 

= a2^t} 
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lienee, we obtain. 


w(x,t) = a(t)x + b(t) 


(3.1.13) 


Ap plying-,' the boundary conditions (3 .1. 12) we havej 
b(t) = a^Ct) 


a(t) 


[a,(t) - a, (t) ]/t 


how, v(x,t) satisiies the partial diiXei'eiitial equation 
9 w(x. t) K_ ( X , t ) _ _ ^w ( x, t ) 


^,^2 at 

whei'e w(x,t) is a lorioiiini function of x and t. 
The initial and boundary conditions becomes 

I.Q . v(x,o) = T(x) - w(x,o) = g 
v(Ojt) = o; v(<^j, t) = 0 

where, w(xjt) is given by, 

w(x,t) = “ [a^ct) - a^(t)]x + a^(t) 

Ther ef ore, 

§w X da^(t) ^ dc:^(t) 

~ gt ~ dt dt" ^ ' dt 


g(x) 


let us use the notation 


dao(t) da, (t) da, (t) 


nx,t) = - 1 (-1^ ^ 

Thus, our equation for v(x,t) becomes 




(3.1.14) 


(3.1.15) 


(3.1.16) 


(3.1.17) 


x-fhere, or= K/pC^, with boundary conditions 


v(o,t) = o 
and initial condition 


(4t) = o 


(3.1.18) 


v(x, o) 


T(x) - xf(x,o) 


(3.1.19) 
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This is an inhomogenoous equation with homogeneous boundary 
conditions. 

Let us seek a solution to this boundary value 
problem in the form 


y 


nnx 


v(x,t) = ^ ^ 

n=l 


(3.1. 20) 


So that the boundary conditions (3.1.18) should be satisfied, 
let us suppose that the function f(x,t), regarded as a function 


of XyCan be expanded in Fourier series. 

f(x,t) = Z-- f^(t) sin ^ 


( 3 .. 1 . 21 ) 


I 


n=l 


f(x,t) sin dx 


nTix 

£ 


where, f (t) = 

^ Q o 

Substituting the series (3.1.20) and (3.1.21) into equation 
(3.1.17), we obtain, 
dT (t) 


fc"-- “ (-^)^ T^(t) ~ f (t)] sin 


nrcx 


-dt 


£ 


n 


n 


= 0 


mza 


Let us put we obtain 


<iT^(t) 

”"dt 


n 




fjt) 


( 3 . 1 . 22 ) 


By using the initial condition for v(x,t) 

c<> 

v(x,o) = Sin = 

n=l ^ 


g(x) 


We obtain initial condition for s 

T„(°) = I f eM Sin ^ 


dx 
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The solution of eq. (3.1.22) becomes 




T^(t) = e [T^(o) + 


J e'^ 


(t-T) 


f^(T)dT] (3.1.23) 


Substituting (3.1.23) in (3.1.20) we have, the complete 
expression for, v(x,t) follows as, 


OO 


V 


(x,t) = ^ / e [Tj^(o) + j e ^ f^(i:)dT]j 

0 


j 

n=l V 


X Sin (^) 

'i ^ 


(3.1.24) 


where, = “ I g(x) sin (^^) dx 


Thus, the complete solution to the problem is the 
combination of both the partial solutions v(x,t) and w(x,t )5 

X e 0 « 


6^ 




^n 


T(x,t) ® 

n=l 


-p^(t-T) 


n 


fj^(T;)di:] ) sin(^^) 


I [a 2 (t) - a^(t)]x + a^(t) (3.1.25) 


Transient numerical Technique ; 

When the transients are severe and the transient 
boundary conditions are of such form that a mathematical solution 
to the problem is not possible. In these cases, the problems 
are best handled by a numerical technique with hand calculator 
or computers depending upon the accuracy desired. 
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In the setup for the problem the second derivatives 
in equation (3.1.1) is approximated by 

^^T(x.t) ^ _J_. 


and 


d^‘ 

2^ 


( Ax) 


- 2 ) 
m+1 ^ V ^ m-1^ 


mn+1 rnll 

3T(x,t) _ \ 

At 


In tliis relation, the superscript designates the time increment 
and subscript corresponds to space increment. Combining the 
relations above, the difference equation equivalent to 
Eq.(3.1.1) is 


rp^ _ p rn^ t 

•^m-H ^ -^m ^m~l ^ 


mn+1 

1. m 


T 


n 

”111 


a 




(5.1.24) 


Thus, if the temperature of' the various nodes are known at 
any particular time, the temperature after an increment of time 
Af niay be calculated by writing an equation like eqn.(3.1.24) 
for each node and obtaining values of T^ . The procedure may 
be repeated to obtain the distribution after any desired number 
of time increments. Equation (3.1*24) can be conveniently 
written as. 


„n+l 

m 


g^t 

Ax^ 


(^m ^m-1^ (5.1.25) 


Ax' 


m 


In this explicit finite difference method an approxi 
matioh associated with difference equation (3.1.25) is that, 
while the temperatures at T^ changes during a small interval 
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to a value the values of t and -, are assxamed 

m m+i m— ± 

to remain constant. 


If the time and distance increnents are chosen such 
that M = = 2, the temperature of node m after 

the time increment is given as the arithmetic average of the 
two adjacent nodal temperatures at the hegining of the time 
increment. The solution of the value of the parameter 
M = ^ o^^ ' t importan.t, as it governs the case with which 

we may proceed to effect the numerical solution. 


Once the distance increment and the value of M are 
established, the time increment is fixed. The smaller values 
of A. X and A t in the iteration process gives better accuracy. 
One important thing to note is that, if M 2 the coefficients 
of equation (3.1.25) becomes negative and a condition is 
generated that the 2nd law of thermodynamics is violated. 

The difference equation given above is useful for 
determining the internal temperature in the core as a function 
of space and time. Once the boundary conditions are specified, 
their finite difference approximations can be made and used 
to start the iteration process. 
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5 • 2 STEADY STATE NATURAL COImTECTION (THEBI40 SYPHON) COOLCTG 

The most important aspect of the present work is 
concerned with the natural convection cooling (thermosyphoning) 
in the horizontal CAKTEU reactor core in the vicinity of each 
fuel rod following the transient, specifically in situations, 
when equilibrium natural circulation in the primary system cannot 
be established and auxilary recirculation pumps become inopera- 
tive. A central question then becomes, whether thermosyphoning 
in the horizontal core plays dominant role in determining 
availability of adequate cooling following the accident, by 
directly extracting heat from the surface of the each fuel rod. 

Problem Formulation s 

A fuel bundle in a CAKTDU reactor core is a circular 
array of circular cylindrical fuel elements placed according to 
a specified array geometry. The array is generally divided into 
triangular lattice cells as shown in Fig, 3. S’; and coolant channe 
is characterized by its' hydraulic diameter, , where, 

A- is the cross sectional area of the lattice cell and P is the 

w 

t 

wetted perimeter. To formulate the problem, we consider indivi- 
dual fuel rods in the lattice cells, which are isolated from each 
other and try to develop an analysis for a single fuel rod. 
Without any ambiguity, the treatment can be extended for all the 
fuel rods assuming that these are all in the same thermal state. 
This seems reasonable because, the stead.y state heat flux 
distribution in the core has been considerably flattened compared 
to the initial sinusoidal distribution. As a result, this gives 
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approximately a constant miform heat flnx throughout the core. 

ViTith the aid of the argument given by Lighthill [7], it 
is recalled that a boundary layer regime will be developed in the 
vicinity of each fuel rod because of the existence of considerably 
hi^ heat flux inside the core. Calculation of steady state 
Grashoff number indicates that the boundary layer regime is 
laminar in nature. Since, the thickness of the thermal and 
hydrodynamic boimdary layers are approximately equal in the fluids 
with Prandtl number close to unity, therefore, no separate 
boundary la^/ers will be considered in the present treatment. The 
region near the fuel rod surface in the boimdary layer regime, 
where velocity is found to be significant. The fluid beyond the 
boimdary layer regime is assumed to be stationary, since the 
buoyant force due to gravity, causes the lighter fluid near the 
fluid to move and it becomes insignificant at a larger distance 
from the surface. The fluid velocity, therefore, be aero at the 
fuel rod surface, increase to a maximum value in the immediate 
neighbourhood, and then decreases asymptotically to aero at the 
outer edge of the boundary layer. The motion of the flow is 
restricted to a region close to the fuel rod surface. For the 
present treatment the following assumptions are made ; 

(1) Mutual influences of the fuel rods on their individual flow 
in the boundary layer regime is disregarded. 

(2) The interaction of the boundary layer of a fuel rod with 
its' neighbouring one is also disregarded. 

(3) The results derived for a single fuel rod is equally 




Fig. Co-ordinate System 


g/5(T-Too) 



Fig. 3-7 Buoyant Force along the 
curved surface 
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applicable to all the rods inside a fuel bmdle according to 
foregoing arguments. 

Usually, Nusselt number and Grashoff nijmber are based 
on the diameter of the cylinder, the significant length in all 
free convection problems with horizontal cylinders. A coordinate 
system suggested by Zoh and Price [22] to treat axially symmetric 
horizontal objects in free convection is well suited for the 
present problem. This helps immensly to construct the buoyant 
force term. The coordinate system is shown in Pig. It is 

located with the origin at the lower end of the cylinder. 

X coordinate varies across the fuel rod surface as shown in the 
figure and y axis is normal to the surface. Here T^ , ^w^^^ 
and T<s^ are temperatures at origin, at a point with coordinate x 
on the fuel rod surface and, ambient fluid (coolant) tonperature 
respectively. 

Isothermal and Uon-isothermal surface conditions 
on the fuel rod ; 

As previously pointed out that, accumulation of volatile 
fission products in the fuel cladding gap and in other void 
regions in the cladding envelope is likely to affect the gap 
conductance and give rise to non~isothermal surface conditions 
on the fuel rod. In such situations, the temperature varies in 
an arbitrary manner with the coordinate x across the fuel rod 
surface. At x = 0, a stagnation point is formed. For the present 
situation, the surface temperature is assumed to vary -in both 
linear and nonlinear manner as two separate cases. The surface 



temperature is assumed to take the following formss 
For linear vari£i,tion } 

(T,, “ To.) = + a(|)] (3.2.1) 

For nonlinear variation ^ 

- tco) = (I„ - T„)[l + a(|)^] (3.2.2) 

0 

V/here, a is an arbitrary constant and the surface temperature 
is alwaj'S higher than tho free stream (ambient) temperature . 

In isothermal surface condition, there is no variation 
of temperature occurs across the surface of the fuel rod, which 
implies that a = o. Hence, it is reasonable to suppose that the 
free convection heat transfer from a non-isothermal surface would 
be significantly different from that of an isothermal one. 

Construction of the buoyant force term s 

The density has been considered as only variable in 
forming the buoyancy force term. Variation of all other ■ 
properties are neglected. The body force is assumed to be 
significant only in the vertically upward direction, then body 
force per unit volume in the vertical direction is gp(T-T^), 
where, g is the gravitational force vector and f. is the thermal 
expansion coefficient. In order to construct the body force 
expression for a curved surface, a tangent is drawn at any 
point B on the fuel rod surface as shown in the Fig. 3.'7. If 
R is the radius of tho fuel rod and x is the coordinate of the., 
point B, then the angle Q subtended at the centre of the circle 




is equal to x/R. The direction of the buoj^'ant force at the 
point B is along tho direction of the tangent drawn at that point 
and its magnitude is gp(T - T^^ ) sin(x/R) . Thus, generalizing 
the buoyancy force at any point on the curved surface of tho fuel 
rod has a direction along the tangent drara at that point and 
magnitude equal to =gp(T - Tg^ ) sin(x/R) , which is a function of 
the coordinate x. 


Governing Differential Equations % 

The two dimensional equations expressing conservation 
of mass, momentum, and energy for steady, laminar flow in a 
boundary layer on a horizontal cylindrical fuel element can be 
expressed, without any loss of generality in Cartesian coordinates 
as follows s 


^ ^ = 0 

dy 


4. -y 

" 3X dj 


5T , ,, dl 
dy 


+ gp(T-a!jsin(|) 

dy 



( 3 . 2 . 3 ) 

( 3 . 2 . 4 ) 

(3.2.5) 


where, u and v are velocity components along x and y directions, 
respectively, -y* is the kinematic viscosity and Pr is the Prandtl 
numb er . 


Viscous dissipation and work against the gravity field 
have been omitted from the energy equation since they are negli- 
gibly small. The appropriate boundary conditions are s 

y=:o u = v= o T = T^(x) 

y— ^60 u — >0 T — Too 





First iDoundary conditions is the no slip condition on 
the surface and second one is the asymptotic boundary condition 
at the edge of the free convection boundary layer. Equation 
(3.2.3) can be automatically satisfied if a stream fimction is 
introduced such that 


u 


di 



(3.2.6) 


The momentum and energy equations may bo converted into dimen- 
sionless forms by introducing suitable similarity parameters as 
shown by Schimdt [23] and Valentine [24] for laminar free 
convection problems. The similarity parameters are given as: 


X = 


X 
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ri 


= [ 


eP(\ - T«) 


.2 


1/4 

] ^ 


R 


= (Gr) 


1/4 Z 
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M(x,ri) = [ 


11 


1/4 


{T„ - T^.)e3 

0 


a 


(Gr) 


1 J1 
^74 V 


0(S,r]) = 


1 

■- 

T - T 


(3.2.7) 


Co 


In the above transf ormations, the first two quantities, i.e. 

X and p represent dimensionless coordinates. M(x,p) and 0(x,p) 
are dimensionless stream function and temperature respectively. ! 

The reasoning behind the selection of this form for the simila- ; 

rity parameters is rather involved and specific for particular 

I 

I 

type of problem. Using these transformation we deduce the following! 
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Substit-uting the above results, the momentum and energy equa- 
tions i.e. equation (5.2.4) and (3.2.5) nia.y be rewritten in the 
following dimensionless form ; 


d¥L 2%. ’M _ 

aian"”^R" 


sin X 


C3.2.8) 


I ^n 2x dx] 


(3.2.9) 


In terms of the new variables, the boundary conditions are 

0 M=g = ^ = 0 (3.2.10) 


h 


o 


0 



and, 0=1 (isothermal case) 

0 = (1 + ax) (linear variation of temperature) 

0 = (1 + ax ) (N on “lin ear ■ variation of temperature) 

n _> oo — ■> 0 0 5. 0 

Solution Procedure i 

Equations (3.2.8) and (3.2,9) together with the boundary 
conditions can be solved by using Blasius-Howarth series approxi- 
mation technique for the dimensionless stream function M and 
dimensionless temperature 0 as discussed by Chiang and Kaye [25]. 
Thus, the series approximations for M and 0 take the following 
forms for non-linear and linear variations of the fuel rod 
surface temperature. For non-linear case ; 

M(x,ti) = X F^(ri) + x-^ F^('n) + .... (3.2.11) 

0(x,T)) = G^^(p) + x^ ^I^b) + .... 

and For linear case ; 

M(x,ri) = X Fq(ti) + x^ F^('n) + .... (3.2,12) 

0(x,q) = ^^('n) + x G 2 (ti) + .... 

Ihe first terms in each of these series expansion is the 
dominating one with respect to the other terms, Hi^er order 
terms are treated as perturbations to the first term. The great 
advantage lies with these series expansions, which convert the 
governing partial differential equation to a set M. 
differential equations. 


I wy l o iiii wi i 

A 6*872 




Sin X in equation (3.2.8) is approximated by the 
first two terms of its power series as follows ; 
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Sin X = X - -rr + . . . . (3.2.13) 

When equations (3.2.11) and sine series approximation (3.2,13) 
are substituted into equation (3.2.8) and (3.2.9) and terms 
arranged in ascending power of x gives. 


dP^ o 

3 - " 0 - 


d% 

dp 
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(3.2.14) 
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dp ' 
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dGj 

o dp 




(3.2.15) 


WM (Y 

How, setting the coefficient of x to zero for each non-negative 
integer a ^3 and neglecting the coefficient of hi^er order 
terms in x, yields the following sets of ordinary differential 
equation 


d3p 

d%' 

dPo 



— ^ + P 
dp^ 
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0,2 

dp 

~ ^ (^p 
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a%, 
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° dp^ 

dP 
^ dp 
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+ Gj^ — 


hr. = 0 

D 0 


(3.2.16) 


(3.2.17) 



0 


1 4 ®1 

d,l2 


dp “ 


dp 

The boundary conditions 

are 


p = 0 

^0 = 

II 

II 

o 

o 

: 1 


= 

II 

o 

li 

: a 

p — ^ oo 

dp 

— ■> 0 

Gq " 

— > 0 


dp 

^0 

Gi- 

0 


(3.2.18) 


Similarly, when equations (3*2.12) together with the sine series 
expansion (3.2.13) are substituted into equations (3.2.8) and 
(3.2.9) and terms arranged in ascending powers of x, the setting 
to zero of the coefficients of x® for each non-negative integer 
a ^2 yields the same sets of ordinary differential equations 
as given above. 

An inspection of equations (3.2,16) to (3.2,18) reveals 
that the present problem involves two parameters, Pr and a. 
However, for isothermal . surf ace condition 'a' becomes zero. For 
non -isothermal conditions, for each combination of Pr and a it 
would be necessary to solve two pairs of simultaneous equations. 
She Prandtl number is considered to be a constant. However, a 
can have different values as it is taken to be an arbitrary 
constant. Hence, it may be desired to solve equations (3*2.16) 
to (3.2,18) for different values of a. Consequently , refering to 
Struble [26], the following transformations are used to transform 
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the original equations and boundary conditions to new sets 
of differential equations and boundary conditions which are 
free of the parameter a. 




Xn (P) 


=Xt(p) 


1 -.. 

F^^Cp) = aX2(p) + X^Cp) &3_(p) = ^2^^^ + X^Cp) 

Substituting these transformations into equations (5»2.16) to 
(5; 2.18), and with little algebra yields the new differential 
equations together with their boundary conditions. 
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dp^ 

dX 
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dq 
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+ 3X, 


d^X 


- + ^3 - 1^1 = ° 
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dY^ 

‘3 d^rJ "" ° 


13 = 0 
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( 3 . 2 . 22 ) 


and hence solutions of these equations are applicable for 
different values of a. 


It is noted that equations (3.2.20) are non-linear, 
while remainder of the equations are linear (assuming that 
solutions of equations (3.2.''20) are known). All of these three 
sets of equations are coupled third order ordinary differential 
equations with asymptotic boundary conditions. Inspection of 
equations (3.2.20) reveals that these are very much similar to 
the boundary value problem for the free convection flow about a 
vertical plate refering to Ostrach [27]. Equations (3.2.21) 
and (3.2.22) involve four dependent variables each, and and 
Y^ are two of these four variables. These asymptotic boundary 
value problems are virtually impossible to solve by an analytical 
method. Extensive numerical analysis of these equations are 
required to solve them fairly accurately by a reasonably good 
numerical technique. The mathematical apparatus required for 
the numerical solution of these equations will be introduced in 



succeeding chapter. 

Velocity, Temperature and Heat transfer i 

Substituting the transformations given by equations 
( 3 . 2 . 19 ) into equations (3-2.11) and (3.2.12), the dimensionless 
stream functions and temperature assumes the following forms 
respectively for non-linear and linear variation of surface 
temperature. 

M(x,ti) = X X^(ti) + x^ (aX2(ri) + X^(ti)) 

(2f(5E,p) = Y^(ti) + x^ (al2(ri) + Y^(ti)) 

M(x,ti) = X X 2 _(ti) + (3X2(11) + X^(ii)) 

0(x,r)) = l3_(ii) + X (ay2(il) + X^(p)) 


( 3 . 2 . 23 ) 


(3.2.24) 


Prom the expression for similarity parameter related to dimen- 
sionless stream function M, given in equation (3.2.7) » we have 



V, . which gives, 

dM(x,r)} 

R an 


Substituting for ^ from the equations (3. 2. 23) and (3.2.24) 
into the above equation, the expression for velocity distribution 
in the boundary layer follows for both non-linear and linear 
cases of thermal condition on the surface s 


u 


u 


/p„\l/2w dX-Y rz dX^ hX-T 

fx —A + x^ fa — - + ) 1 

R 1-^ dp dp ^ dp 


R '•dp dp dp 


( 3 . 2 . 25 ) 
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Now, equating Fourier law of heat flow on the surface in terms 
of 0 and heat flux, q/A given as. 


q/A = - K{I^ - T^,) 




■o ® "''1,1=0 

and Newton's law of cooling for the heat converted away from the 
surface given as, 


q/A = h (T, 


w. 


^ Cks ^ 


One finds the follox^ing equation 
h = 


„ j 

3’i|,=o 


The dimensionless heat transfer i.e., Nusselt number may be 
given as. 


= ^ = 0 (3.2.26) 

where, NUjj based upon the radius (i.e. diameter) of the fuel rod. 

Substituting for I ^ from giequations (3.2.23) and (3.2.24) 

cfT] f ri=u 

into the above equation, the expression for Nusselt number follows 
for both non-linear and linear cases of thermal condition on the 
surf ac e 


Nup = - Gr^/"^ 


. .-V 

dt) 


dY.(O) dY,(0) 


dp 


dp 


[ 


+ X (a 


dY.,(0) dY,(0) 


dp 


dp 


)] 


■)] 


(3-2.27) 


She equations representing, velocity and Nusselt niamber for 
non-linear and linear cases are similar except for the power of x 
involved in those expressions. The velocity distribution in the 
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■boundary layer accounts for the thermosyphoning effects, in 
the form of eddies in the vicinity of each fuel rod. In order 
to compute the velocity and Nusselt number, the derivatives 
in the corresponding expressions must "be evaluated which are 
functions of t) , the measure of the dimensionless distance from 
the surface of the fuel rod. 
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CHAPTER r\?- 

RIMER IC All ANALYSIS OP THE SYSTElvIS 
OE BOUNDAJIY LAYER DIEEERENTIAL EQUATIONS 

4.1 INTRODUCTION 

The ntimerical integration of the ordinary differential 
equations, obtained from Blasius-Hawarth series with appropriate 
transformations, involves the satisfaction of asymptotic boimd- 
ary conditions. That is, some boundary conditions are specified 
at the initial point or wall and other are specified as limits 
that must be approached at relatively large values of the 
independent variable, corresponding to the edge of the bomdary 
layer. In order to integrate the differential equations 
nmerically from the wall to the edge of the boundary layer, it 
is necessary to specify as many additional conditions at the 
wall as there are conditions to be satisfied at the edge of the 
boundary layer. All methods that have been used to find the 
proper initial conditions rely on the fact that, for relatively 
large values of the independently variable, the integrals of 
the differential equations depend on the initial conditions. 
Among the broad categories of numerical methods, two widely 
applicable techniques are namely : 

(i) Initial value method, Eox [28]. 

(ii) Quasilinearization method, Radbill [29]. 
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Initial value method, that has "been used to find proper 
initial conditions, is that of obtaining integrals of the diffe- 
rential equations with guessed initial conditions, that is an 
attempt is made to integrate the differential equations to the 
edge of the boimdary layer. If this can he accomplished, 
corrections are made for the initial guesses hy suitable itera- 
tive method and the process is repeated untill desired accuracy 
is achieved. In order to carry out these iterations, additional 
differential equations have to be integrated. The additional 
differential equations are obtained by differentiating the terms 
of the original differential equations with respect to initial 
conditions. The integrals of these additional differential 
equations, referred to as perturbation equations, give the rate 
of change of the integrals of the original differential equations 
with respect to the initial conditions. 

Quasilinearization is similar in principle to the initial 
value method, except that the original differential equations, 
if non-linear, are linearized. The resulting linear differential 
equations for the current approximation are inhomogeneous, since 
they also contain the previous approximation as members. In 
application, the two methods are different in that, in the 
quasilinearization method, the complementary solutions are used 
to obtain solutions to the inhomogeneous linear differential 
equations, whereas in initial value method, the integrals of 
the perturbation equ^'tions are used to adjust the starting 
values at the initial point. 
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A disadvantage associated with quasilinearization 
method is that, the solution of the inhomogeneous equations 
obtained by combining complementary solution and a particuilar 
integral is usually not well determined except near the origin 
or initial point as pointed out by Hartree [50], Thus initial 
value method has an edge over the latter, as solutions of pertur- 
bation equations can be used to determine the proper initial 
conditions closely. Since, boundary conditions specified as a 
limit that must be approached at higher values of the independent 
variable, the integration should be stopped at a value of the 
independent variable, so that various functions are approaching 
their asymptotic and physically acceptable values. This 
value of the indopondont variable may bo callod approximately 
t ho edge of the boundary ..layor. This value of independent 
variable is obtained by suitable guess and trial effort. 

In an effort to adapt the initial value method to the 
solution of differential equations with asymptotic boundary 
condition, an iterative method developed by Uachtsheim and 
Swigert [31] to correct the initial conditions, is applied. 

This method of solution is capable of satisfying boundary layer 
correctly. This is accomplished by choosing the additional 
initial conditions so that the mean square error between the 
computed variable and asymptotic variable is minimized. This 
method of solution is based on least square error criterion and 
the edge of the boundary layer is approached in steps. 
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4.2 DESQRIPTIOF OF TFFl METHOD 

A description and application of the initial value 
method can best be given, refering to any one of the set of 
differential equations as an example. However, the general 
iterative formulas derived by using first set of non-linear 
equation (3.2.20) are equally applicable to all the sets of 
equations. 


Using 'primes', to denote differentiation with respect 
to independent variable p , the measure of the dimensionless 
distance from the wall, we have; 


dX. 

dp 


d^X. 


= X 


1 ' ,2 
dp 


d^X-, 

■y — IT ^ 

A-j ^ •:? — A-j » 

dn^ 


Nowj, set of equations (3.2.20) can be rewritten as 


VHI 


t — 


- (x^ x" - X^^ + I^) 

- Pr Xi 


(4.2.1) 


with boundary conditions; 

n = 0, Xj^Ct)) = = 0, = 1 

XiCn) = x^(,i) — 0 

The dependent variables, X^ a^d related to dimensionless 
stream function and temperature respectively. 

In practice, the asymptotic boundary conditions is 
replaced by conditions that X^ = 0 and = 0, to a sufficient 
degree of accuracy at p = where is value of 
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the independent variable. IThe boundary value problem of this 
coupled differential equation is equivalent to the problem oi 
finding values of X^(0) and y^(0) for which the boundary conditi 
ons at the edge of the boundary layer is satisfied, i.e. it is 
desired to find solutions of equations 


< edge = 0 


^1 edge tXi(O), y;(0)] 

) 


= 0 


( 4 . 2 . 2 ) 

(4.2.3) 


where, = x' (p 


edge “ "^1 


edge^ * 


1 ''''edge^ 

Since there are two asymptotic boundary conditions to satisfy, 
two additional initial conditions X^(0) and 12^(0) at the wall 
will have to be adjusted, that is, values of these are sought 
that will satisfy simultaneously the above equations. Tbe 

functions X^ gdge^^l^^^ ’ y][_(0)] and edgeC<(°)- ^1(05 ;i ^ 
this problem are not, in general, explicit; those will be 
expressed as functions of 21^(0) and throu^ integration 

of equation (4.2.1). 

Using notation^ x = X^(0) and y = Y^(0), we observed 


that, small changes /\ x and X\y in the initial conditions change 


li_ 4ix + 


and /^y 


X^ and Y^^ "by the amounts -r ewiu. 

respectively. The necessary correction to the first approxi- 
mation comes from the solution of the linear equations 


, ax' d'^\ 


= 0 


d^x 


3^1 ^ 


= 0 


(4.2, 4) 
(4.2.5) 


at p = p 


edge" 


* TbA nntat.iOYis .hRrfi should not bo nonfuSAd with the? nonrd inatss- 
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The initial value method again demands that all higher 


derivatives of the independent variables vanish at the outer 


boundary i.e. at p = 'Hedge* Hence, satisfaction of both pair 


of boundary conditions i.e. = 0, = 0 and = 0, = 0 

at p = 'HgjjgQ ensures that Xj_ and approach zero asymptotically 
as p approaches infinity. Thus, in order to satisfy asymptotic 


conditions at a finite value of p, both pair of above conditions 


must be imposed. Hence, in order to satisfy, asymptotic 


boundary conditions, the preceeding equations (4.2.2) and 
(4.2,3) must be supplemented by 


edge = 0 

edge K <0). (0)] = 0 


( 4 . 2 . 6 ) 


( 4 . 2 . 7 ) 


Again for small changes Ax and Ay in the initial conditionSs 
the necessary correction to the first approximation, as 
previously given, comes from the solution of linear equations 

« , ax" . , . 


« d^'i ax" 


( 4 . 2 , 8 ) 


• ^^1 A ^^1 


( 4 . 2 . 9 ) 


at p = p 


edge* 


The solution of both pairs of equations i.e. (4. 2 .^ 3 , 
(4.2.5) and (4.2.8), (4.2.9) can be performed, provided that 
the partial derivatives involves. in those equations with respect 
to X and y can be evaluated at p = ®He partial deriva- 

tives can be evaluated by forming and integrating perturbation 
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I H t 

equations for the functions Xj^, X^ and Y-j^, . These 

equations are obtained by differentiating the terms in the 
original equations (4.2.1) approximately with respect to x and 
y respectively. 


The perturbation diff er«ntial equation for x deriva- 
tives are ; 


3x' 


(^X 


<5Y. 
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, 8x"v , 3J\ 

<-7r h + * ^iw 


32i . 3^1, 

Pr (,^ + X3, g-^) 


(4.2.10) 


with the initial conditions 

n. !5= Ei ?Ii 0 

^ ~ ^x ~ 3x“ ~ dx “ 0X “ 9x 
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The perturbation differential equations for the y derivatives 


are 


ax 

w 


1 


j)x „ ax’’ , ax^ 


§nd 


dY 


«) 




ay 


3Xi , aY^, 

^ay ^1 % ay ' 


?y 


(4.2.11) 


with initial conditions 

ax^ ax^ ^ ^ _ aY^ 

ay “ ay~ ^ ~ 3y ~ ay “ 


Althou^, the equations for the y derivatives are the same as 
the equations for the x derivatives, but the integrals . • 

of these equations will differ since the initial conditions are 
different. It is also noted that there are as many system of 
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perturbation equations to integrate as there are asymptotic 
boundary conditions to meet. 

Uow, corrections to the initial conditions i£i>x and Ay 
so chosen that these four equations (4.2.2), (4.2,3), (4.2.8) 
and (4.2.9) are satisfied at q = This is , in general, 

impossible since there are four equations and only two adjusta- 
ble parameters. However, a satisfactory solution that is 
consistent with the idea that the boundary conditions can't be 
satisfied exactly at a finite value of q , is to seek least 
Square solutions of those four equations. 

4.3 HACHTSHEIM-SWIGERT ITERATION SCHEME 

The method developed by Hachtsheim and Swigert [31] 
can be applied to these four equations to solve for Ax and Ay 
basing on least square error criteria. According to the method, 
Ax and Ay should be found such that, the sum of the squares 
of the error in these four equations must be minimum. Hence, 
c'oasiderang discrepancies 4n these equations 

respectively (deviations from their asymptotic values), where, 

I 

^1 ^ ^x ^y ~ 

„ dx" 

^1 = 

, axi 

It is equivalent of saying that, the original equations (4.2.2) 


(4.3.1) 
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and (4.2,3) along with the snpplenented equations ( 4 . 2 , 6 ) and 
(4.2.7) are associated with these descripancies 6 ^, b 2 ^ 6 ^ and 
6 ^ respectively. We have s 


edge ( 0 )> ( 0 )] = 

■^1 edge H (°>] = 

edge ( 0 )- = 

edge (0)3 = 



(4.3.2) 


Solution of set of equations (4.3.1) in least square sense 
requires detemination of nininun values of the error, 

E = 6 ^ + 62 + 6 ^ + 6 ^ 

So differentiating E with respect to x and y respectively and 
equating to zero i.e. we evaluate; For x derivatives, 


il = -1 

3x dx 


(&l + 6^ + 6^ + 6^) 


2 2 

0, Substituting for ^ 2 * 


0 o 

6 ^ and 6 ^ fron the set of equations ( 4 . 3 . 1 ) » we have after 
rearranging terns; 

r(!h)2 ^ (Ei)2 ^ ^ ^ ^ + 

L(3x > (jx ) + (gx ) + (jx ' J + (ax ay dy 


n If 


dF‘ ^ fF 


, asq 31 1 „ 

•*-. . -rr ^ JU . V 4- 


■] /sy = - (X^ -jA + Tj + X 


I, =!xr(") 


Sinilarly, for y derivative, equating to zero, and 
rearranging terns, 


"1 ax^ ^1 ^x 

(4.3.3) 
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rEi + El El + til Eim^ ^ rfEiis ^ (£1,2 

‘-^x §y ^y. 5x 3y ^x §y ‘■^0y ^ 


+ (• 


3 X 1 2 ?Y. 2 

i)^+(^)^]Ay = 


di 


- (z 


az, 3y „ 3z 
1 3y ■*■ ^1 3y ^1 m 


5Y. 


"1 3y 
(4.3.4) 


Equations (4.3«3) and ( 4 . 3 , 4 ) can be written as a system of 
linear equations in two unknowns in a simplified form. 


= ®1 

■" = ®2 


(4.3.5) 


and in matrix notation. 





^yy 


y 


/b. 


\®2 


where, 

■^x 


az 


If 


1,2 . El, 2 . El, 2 . El, 2 


^ + < 0 r> + 


3X, ax. 3 I-, 31. 9X., 3X-, 3Y-, 31. 

^y iF” ^ W ^ Wr '*' P"" 


0 X -1 ^ 'hY Q TO V j. n ^ 

V = ^ ‘ 3 r> + ^ 5 r> 


■ 3X. 


9y. 


B- 


, ax. 




1 3 x 


i »i ax. I c/Y. 
- + X. ^ + Yt 


^1 32 : 
„ ai.' 


1 3x 

, ^y- 


r ax. 3Y. „ 1/.41.. , t/A. 

sa = - <^i aT + ""i ^ ^1 3^ ^ ""i 

Solving equations (4.3.5) for the two unknowns Ax and Ay, we 
have. 



oy 


and 


Ax = ~ 

■*^xx ■*yy"^y 


A y 


^2 ~ ^1 
■^xx \y"\y 


(4.5.6) 


(4.3.7) 


This calculation yields values of A x and Ay corresponding 
minimum of the sum 6^ + 6^ + ^ ^ 4 * magnitude of error E, 

evaluated for a given range of independent variable approached 
in steps, gives an indication how unsatisfactory the asymptotic 
boundary conditions are. The error teimi can also be written as 

I P p tip t p 

E = (4.3.8) 


The minimum of E with respect to x and y, say corresponds 

to no deviations in the required initial conditions, that permits 
the asymptotic boundary conditions to be approximately satisfied 
at a finite value of the independent variable. The edge of the 
boundary layer, taken to be some guessed value of r) consistent 
with our assumptions, for which is less than some preassigned 

small value for this range of integration, gives the desired 
solutions. 

These modifications to the initial value method come 
from the corrections to the initial conditions given by equations 
(4.3.6) and (4.3.7) in a least square sense. These are genera- 
lized correction formulas and can be applied to any system of 
coupled third order differential equations, if the subscript *1' 
associated with X and Y:feignored. However, the main difference 
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lies with the evaluation of corresponding partial derivatives 
involved in correction equations, which would he different for 
different set of equations and also correspondingly their 
perturbation equations. 

Application to systems of equations (3.2,21) and (3.2.22) ; 

In primed notations as stated before, these equations 
can be rewritten ass 



In the above equations, dependent variables each. These 
equations can only be solved by knowing the solution, namely, 
the values of corrected initial conditions X 2 ^( 0 ) and of 
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equation (2.4.1). Method of solution of these equations is 
exactly same as given in the foregoing analysiSj once the values 
of the initial conditions are evaluated. The only difference 
lies with the evaluation of partial derivatives involved impli- 
citly in the expressions (4.3.6) and (4.3.7), which would he 
different for different system of coupled equations depending 
upon the integrals of corresponding perturbation equations. 

The perturbation differential equations for the 
coupled equations (4.3*9) are obtained by differentiating 
original equations with respect to the initial conditions x and 
y, respectively. Por x derivative. 






'3 X 


0T, 

< 

2fx 


dx, „ mp , $Xp 

- (^^ X 2 + ^ 4X^ 


^X. , 


It 

ax. 


ft 

3X2 

ax 2 


X2 + 

3 X^ 

Qx 

“ ax ^ 


axo 

1 

ax^ 


= - Pr(X 3 ^ 

ar -^^2 

dx 

2X1 

ax{ 
" ^ ax 


3X2 

dx{ 

0 x 

ax , 
+ 3 


(4.3.11) 


With initial conditions: 


•q=Oi 


TI-OS 


axp 21X2 


Qx 

^x. 


ax' 


ax ax 


d^2 

0x 

m 

§x 


dY]_ 


= 0 , 


0 , 


dx, 

ax' 


1 _ 


For y derivatives, 
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dJ. 


I?f 


, ax 


ay 


2 


^2 + ^1 f7” “ ay^ “ 4 ^ x^ 


ax: 


3^Y, 

9y 


■1 n ax-p ax p 

^ ay 2 + ^^1 3y 5x ^ 


, 07 P , 9X-, aY.’ 

- 5r + ^2 ^ + 3 x 2 + 3 ^ Y^) 


( 4 . 3 . 12 ) 


_1 , , ^ 


With, initial conditions^ 

r)=0; 

r|=0; 


ro 

ax’ 

dll 

-?I2 - n 


dj 


"W~ 

ay ’ 

9y 

^Ii 


_ ?Ei 

3Y, 

= = °> 

ax; 

^y 

ay 

ay 

ay 


= 1 


= 1 


Similarly, the perturbation differential equations for the 
system of coupled differential equations (4.3.10) are as 
follows; 

For X derivative. 


ax 


!SI 


ax 


2 -. _ 


‘ax “3 

ax. 


^ , f 1 (J^ry t 

( 5 ^ ^3 ■*■ 3x " 


t ^ 

ax 


- ^^""3 


■1 . ax, 3X, 3Y, 

+ 3^X3 + 3Xi^+3-^-i|ji) 




, dK . ' 3 x 3 
^ ^3 ^ 


ax, , BY, 

_ p — w4 Y «. PY —2. 

^ ^ 3 1 


BYi ax, 

+ 3 X 5 + 3 Yj^) 


(4.3.13) 


With initial conditions 
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ri=0: 


ri=0'i: 


dX 2»X_ SY 

3x ^x ax 

ax^ 'dY^ 


T^x ax ax 


= 0 , 


ax 


ax 

ax 


= 1 


1=1 


For y derivatives; 


ax 


f !f 


ax. 


ax, 


ay 


_ W-“-T tl U-^'-'Z 

2— _ _/■ i Y a. Y 2. _ 


“ “^ay"’ ^3 &y 


, ax, axn , 
«1 aT - ^ jT ^• 


ax, „ ar, >T. 

+ 3 x, + 3Xi + 3-^ - 5 


u 

ay 


aY, , ax, ‘ax. 

-Pr(Xi ^ ^ - 2 ^ 

ail ax. 


Y_ - 


, ax, 

PY 2. 

^1 ay 


(4.3.14) 


With, initial conditions; 


■n=C: 


n=0; 


ay “ay “ ay “ ay 


ax^ dx, 

3ij“ = 


ax 


ail 

ay 


= 0, 


= 0, 


111 

0y 

ail 

ay 


= 1 . 


As stated previously, the partial derivatives with 
respect to x and y, that appears in the generalized correction 
equations ( 4 . 3 . 6 ) and ( 4 . 3 . 7 ) can he evaluated for these system 
of equations hy integrating corresponding perturbation equation 
for the functions Xg? X 2 » X 2 » I 2 ^ 3 * ^ 3 » ^3 

respectively. 
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4.4 COMPUTATIONAL PROCEDURE 

Usually, a boundary value problon is recast ed as an 
initial value problem by the application of initial value method. 
Hence, the boundary layer and corresponding perturbation diffe- 
rential equations were rewritten as equivalent systems of first 
order differential equations. With proper initial guesses these 
were integrated by fourth order Runge-Kutta method for first 
three values and subsequently, Adams-Moulton predictor-corrector 
method was applied for one correction per step and a fixed 
increment. The corrections to the initial conditions Ax and Ay , 
which involves integrals of perturbation equations were deter- 
mined. The process was repeated until the relative changes 
in the correction terms ^x and Ay are less than small pre- 
assigned values. Least square error was found out and tested 
with a small pre-assigned value. Previous guesses were correct ed 
and the whole process was repeated ■until the desired accuracy 
was achieved, i.e. the least square error was found less than 
a small preassigned value. 

The set of equations (4.2.1) and its’ associated 
perturbation differentials equations (4.2.10) and (4.2.11), 
which are all third order in t) , can be expressed as an equiva- 
lent systsn of fifteen first order differential equations with 
five equations coming from each set. Thirteen initial conditions 
in total are known, and two more conditions are to be guessed as 
initial estimates. These equations subject to the given numeri- 
cal method, lends to the solution within desired accuracy. 
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The two other sets of coupled equations of the 
boundary layer respectively (4.3.9) and (4*3»10), with their 
corresponding perturbation differential equations (4.3*11) > 
(4.3.12), (4.3.13) and (4*3.14) can be expressed as an equiva- 
lent system of twentyfour first order differential equation 
for each set. This is because each of the set of equations 
constitutes four dependent variables. Because of the involve- 
ment of the dependent variables and in these equations, 

the number of initial conditions required to start with, are 
twenty eight and two conditions are guessed as before for each 
set. These initial conditions comprise all the given conditions 
required for the solution of equations (^•;^,1), besides their 
own initial conditions. The solutions of equations (4.2.1) 
achieved by the method just outlined in the foregoing paragraphs, 

tt t 

one had in particular good approximations for X^(0) and Y^(0) . 
These values are then used as a complete information of the 
initial conditions required for the subsequent solutions of 
equations (4.2.10) and (4.2.11). The computational technique 
is exactly the same for these equations. 

The method described above W;„s programmed in FORTRAM 17 
on the DEC -10 Computer system, A general flow chart (Fig. 5.1) 
is presented for the main program adapted for all the system of 
equations. The flow chart of integration subroutine is presented 
in the Fig. 5.2. The FORTRAKT listings of the method are enclosed 
in the appendix D. The integration formulas, used in the inte- 
gration subroutine are also given in Appendix B. 
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4. 5 RESUITS 




and 


With trial initial guesses 0.60 and -0.50 for 5- 

dY^(o) 6.r]‘^ 


dr) 


respectively, the set of equations (4.2.1) were 


integrated to a value of t) equal to 5.0 with initial step 

size of 0.05 for first four steps and subsequently raised to 
0.1 for further steps. The preassigned values to compare least 
square error and increments in the initial estimates in the 
iteration process was kept 0.005 and 0.005 respectively. The 


final results are quoted in Table II. Erom this solution, we 

d2X^(o) 


have in particular better approximations for 


0.796700 


and 


dY^(o) 


dp 


dp 


-0.438100. The least square error was found out 


to be 0 . 000972 . 


These values are utilized further to integrate the other 
two sets of equations (4.3.9) and (4.3.10) with the initial step 
size of 0,05 for first four steps and subsequently raised to 0.1 
for further steps, for the same range of the value of p 




edge 


u pv dY 2( 0 ) 

i.e. 5-0. The initial values ^ — and — ^ — , corresponding 

dp ^ 

to the solution of set of equations (4.3*9) a,re 0.004316 and 

dfe^(o) 

- 0.054234 respectively. The initial values of — and 

dY,(o) dp 

— corresponding to the solution of set of equations(4.3.10) 

are 0.008510 and 0.001510 respectively. 


Solutions of equations (4.2.1) give approximately the 
velocity and temperature profile and their corresponding gra- 
dients in the adjacent boundary layer. These are plotted in 
figures 4.1 through (4.4) respectively. The perturbations to 
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the velocity and temperatures comes from the solutions of sets 
of equations (4.3.9) and (4.3.10) are plotted in figures (4.5) 
through (4.8). 


The set of equation (3.2.27) a.re rewritten with the 
value of Grashoff number substituted from the data given in 
Appendix III as follcrtvs; 


P 1 /. dY-,(o) 5 dY^(o) dY-(o) 

BUj, = -(0. 4262X10®) V4 + ^2;^ ^ -fc—) ] 


for non-linear variation of the temperatures 


dp 


dp 


(4.5.1) 


Nu 


‘D 


p T dY-,(o) dY 5 (o) dY,(o) 

= -(0.4262x 10®)1/^[— ip- + 5(a — |j— + —I—)] (4.5.2) 


for linear variation of temperatures 


NUj^ 


P T / , dY-i (o) dY,(o) 

= -(0.4262X10®)^/^ [-4r-+ ^ -4:—] 


dp 


dp 


(4.5.3) 


for isothermal case, when a = o. 


Now, the dimensionless heat transfer, i.e. Nusselt number 
can be computed from the above equations utilizing the solution 
of boundary layer equations with respect to the dimensionless 
distance. Three cases of .nonr-lsotbDrnaX.-conditiJDns .aye 
considered with three different values of a taken as 0.1, 0,5 
and 2.5 respectively. The results are tabulated in table III 
and (IV) and corresponding plots are displayed in figures (4.9) 
throu^ (4.12). 



Numerical Resvilts of the solution of system of equations (4.2.1) 
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Table III 

Nusselt Number for isothermal and non- isothermal 
surface conditions on the fuel rod 


(a = 0,0, a = 0.1) 


Dimensionless 
distance from 
the origin 

; X 

Isothermal 
surf ac e 
condition 

Si “ 0*0 

' Nusselt number for non-iso thermal] 
surface conditions a = 0,1 J 

Temperature 
vation of the 
f orm 0 = l+aic 
(linear) 

Temperature 
vat ion of the 
form 0 = l+ax^ 
(non-linear) 

0.0 

35.399 

35.399 

35.399 

0.1 

35.393 

35.609 

35.415 

0.2 

35.387 

35.822 

35.479 

0.3 

35.380 

36.043 

35.587 

0.4 

35.371 

36. 256 

35.739 

0.5 

35.365 

36.472 

35.933 

0.6 

35.357 

36.688 

36.171 

0.7 

35.349 

36.905 

36.451 

0.8 

35.341 

37.119 

36.774 

0,9 

35.328 

37.336 

37.142 

1.0 

35.315 

37.551 

37.551 

1.1 

35.301 

37.788 


1.2 

35.286 

37.982 

38.512 

1.3 

35.269 

38.198 

3^.151 

1.4 

35.249 

38. 414 

39.881 

' 1.5 ' 

35.223 

38.629 

40.491 
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Table IV 

Nusselt number for non-isothermal 
surface conditions on the fuel rod 


(a = 0. 5 , a = 2 , 5 ) 


< 

[Dimensi onless 
distance from 
;origin 

X i 

Nusselt Number 

a = 0. 1 

5 

1 a = 2.5 

Temp orature 
variation of 
tho form_ 

^ = (1+ax) 

(linear) 

Temperature 
variation 
of the form 
9f=( l+ax^) 

(Non-linear) 

Temperature 
vari ation 
of the fo^ 
^ = (1+ax) 

(linear) 

! Temperature 
variation 
of the_form 
0=(l+ax^) 

(Non-linear 

0.0 

35.399 

35.399 

35.399 

35.399 

0.1 

35.837 

35.437 

36.517 

35.505 

0.2 

36. 281 

35.571 

37.649 

36.845 

0.3 

36.724 

35.793 

38.777 

36.409 

0.4 

37.168 

36.103 

39.905 

37.198 

0.5 

37.611 

36.502 

41.033 

38. 213 

0.6 

38.055 

36.990 

42.161 

39.452 

0.7 

38. 498 

37.567 

43.288 

40.120 

0.8 

39.942 

38.232 

44.416 

42.612 

0.9 

39.385 

38. 986 

45.544 

' 44.529 

1.0 1 

39.829 

39.829 

46.672 

46.672 

1.1 

40. 273 

40.760 

47.800 

49.040 

1.2 

40.716 

41.781 

48. 928 

51.634 

1.3 

41.160 

42,889 

l 

50.055 

54.454 

1.4 

41. 603 

1 44.081 

i 

51.183 

57.499 

1.5 

42.047 

1 45.328 

i 

52.311 

60,770 



Fig. 4-1 Veloc 




84 



-I — JDeijs ssaiuoisuauiiQ 




Fig. 4-2 Distribution of shear in the boundary layer 
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Fig. 4-6 Perturbations to the dimensiontess temperatu 
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Fig. 4-7 Perturbations to the dimensionless 
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Fig. 4-10 Nop thermal surface condition 
on e fuel rod for a = 0-1 
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Non -isothermal surface condition 
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Fig. 4 -11 
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i‘g.4'12 Nonisothermal surface condition 
on the fuel rod for a =2-5 



PI ^USSIOH OF RTilSTTT.'P .q MD CONCLUSION 


The preponderance of the preceeding chapters are 
devoted to study the problem of transient flow and temperature 
four associated with loss— of— flow accident in a CAUDU 
hedctor core and principally following the transient, the 
tiu.'nnosyph on (Natural convection) cooling effects are assessed 
in situations, when there is a variation of temperature on the 
luul rod surface. 

The results loading to a flow coastdown curve, 

I'ig. 2 , 2 , represents situations when pump flywheel inertia is 
large. Consequently, the pump will continue to circulate the 
fluid for a substantial length of time and the flow coastdown 
will be relatively slow, Following the total loss of p\;mping 
power, the driving pressure is found incapable of establishing 
natural circulation in the primary circuit. 

The t (smperature distribution in the core under the 
transient condition, induced by loss-of-flow accident, 
assuming an initial sinusoidal heat flux distribution can be 
known from the solution of the "transient heat conduction 
equation as cited In the equation (5.1.25). The steady state 
temperature distribution of the core follows from the same 
equation, when time is considered w«ry large. 
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The natural convection thermosyphon results leading 
to velocity and temperature profiles show an asymptotic decrease 
as the edge of the boundary layer is approached. The thermo- 
syphoning elfects are attributed to the velocity distribution in 
the botmdary layer. Nusselt number is calculated for natural 
convection around horizontal cylindrical objects by Hyman 
correlation as givoi in equation (1.3.10). It gives, lluj^=37.38. 

The present heat transfer results as given in Table^J^JE 
Tur iecthormal and non-isothermal surface conditions are in 
close agreement with this value. • Furthermore, the present 
results are approximately one fifteenth of forced, convection 
Nusselt number (for coolant velocity 27.2 ft/sec) as shown in 
the text by K. Sri Ram [32]. The variation of surface tempera- 
ture of the fuel rod in the manner (1 + ax^) gives a sharp rise 
in the Nusselt number compared to the variation in the form 
(1 + ax) for tho values of x '^1, as shown in Fig, 4.11- 

5.2 CONCLUSION 

On the basis of the analysis given in the foregoing 
chapters, the following principal conclusions can be drawn: 

(1) The hydraulic modelling to determine the transient 
flow b«*«wiour asBooiated wltb the loss-of-flow accidents for 
coolant puopa dealgnod with suffioientu large flywheels shows 
that, the flow oosstdown will he slew. However, if the reactor 
trip is iomediate, ths after heat from the decay of fission 

products can be transported out of the core for a substantial 

length of tiao* 
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(2) The steady state temperature distrihution in the 
core is attained approximately after 2 minutes. 

(3) Since, the natural circulation can not he established 
in the primary flow circuit, thermosyphoning in the vicinity 
of each fuel rod provides immediate cooling. To investigate 
this phenomena, following lighthills’ [7] approach, the bound- 
ary layer differential equations for free convection flows 
around horizontal cylindrical fuel elements are formulated and 
solved by numerical methods. For non-isothermal conditions, 
the surface temperature of the fuel rod is assumed to be varied 
in two different manners, i.e. of the forms (1 + ax) and 

(1 + ax^) . Suitable transformations are used so that the 
solutions will bo independent of parameter a and hence, the 
method will be applicable for all reasonable values of a. 
Purthurmoro, it is found that the surface temperature variation 
has a significant effect on heat transfer. 

9. '5 SUGGEST IONS FOR fWTimR . mi 

The present work can be extended in many ways as 
quoted below s 

(a) The hydraulic transient analysis can be extended for very 
fast transients, where acoustical phenomena becomes impor 
tant or to situations where transient boiling and rapid 
depressurization of the systems cuuse rapid time changes 
in the fluid density to occur. These can be accomplished 
by formulating appropriate mass, mcanentum and energy balance 
eqp^tions to such situations. 
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conceive of the pure convection process being relevant 
for the random nature of the flow in the boundary layer 
regime and possibly, neighbouring boundary layer jjater- 
actions becomes prevalent. This may also be investigated. 
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■APPEND TY A 


Some aspects of solution of equation (2.2.10) 


EcjUdliion (2.2.10) can "be wnit'ten as 

1 


H + y2 


(1 + t/tp) 


(A.l) 


subject to initial condition, 

S' = 1 at t = 0, i.e. 

Equ.’ition (A.l) may bo reduced to the form 


y(0) = 1, where, y = ^ 


* 2 . = 0 


(A. 2) 


by i'ho substitution y = zt^ , where f(t) = 

P 

Equ?-i-tion (A. 2) is a special case of Riccati equation 

Po<^)t|t 2^3 + P3_(t)Z + P2(t) = 0 (A.3) 

whore, ■^o “ ' ^1 * ^ ^2 “ 

A theorem, Brand [19] states that the Riccati equation (A.3) 
arid the associated linear equation 


p^(t) + p^(t) || + P2(t)u = 0 (A.4) 

dt 


have corresponding solutions, which satisfy the equations 
u(t) » exp J^z(t)dt, Z(t) = ^ 

By tho above theorem, associated linear equation of (A. 2) is 


found out to bo 


d^u 


dt*^ 


f{t)u 


0 


(A.5) 


with transfonnod initial condition 


1 , . du( t l j - i- 
uTtT at ft*0 "tg. 

Equation (A.5) is of the form of Schrodingr equation. There :k& 
is no standard method of solution of this, equation and usually 
approximation thooiy is used for analytical solution. 
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APPBEDIX B 

ffllSGRAlIOIT F OBMTTT.A.q 

Runge-Kutta 4th Order Pormulas s 

A pth order differential equations caa be expressed as 
a system of p 1st order differential equations, i.e. 


dxP 


ffx.y, ai dJ-V, 

dx dx^^ ^ 


(B.l) 


aP-\ 


It WC toko, y, = y, y^ . g ^ 

P dxP ^ 


the oystem of equations, i. 

0 • 9 

dyi 

dx ~ ^1 •••• 

•’ yp> 

dy^ 

3x""'' ^2 * • ♦ • • 

yp) 

A 

3x “ (^ty> ••••■ 

•’ 

mmmlh 

Voctorially, ^ *= f(x,Y) 

where. 


where, Y, f are column vectors. 
For a couplod equation of the form. 


dy 

dx 


f{3E,y,z) and || = g(x,y, 2 :) 


the 4th ordor Runge-Kutta method is given as 

^n+1 * i ^2 + 2K5 + ^^4) 


^+1 “ A 6 ^ ^1 2^2 + 2^^ + A 


(B.2) 


whoro. 


h = “<*0 . Jh' 
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= he y^, 

K5 =a tlf (x + y + — i 7 4. ') 

2 ' n 2^ ^ 2 ^ ^ 2 ^ 

K.> = hf (x + ^, V + ~Z 7 , \ 

$ ^ n 2» ^n ^ 2 » + 2 ^ 

& = hg + |. Sn + ?, Z^ + 4 ) 

^3 = hg (x^ + |, y^ + z^ + f ) 

= hf (x^ + h, y^ + K^, ^ ^ 

,'4 = hg (Xj^ + h, y^ + K3, Zjj + % ) 

wheru, h io the step size. Three values are first determined 
from the above method to use them in Adam-Moulton predictor 
corrector for hotter accuracy, as corrector 'an implicit itera- 
tion technique is used for next predictor value. 

Adams*-I''ioulton Predictor Corrector Schemes 
Predictor : y^°) = + |j[55f - 59f 

=V,_l) + 37f(V2'yn-2>V2^ " ^n-3’ 



z„ + 
n 


IjESSgCx^.y^.z^) - 59f(=Si-l>yn-l'^-l) + 37g(x^_2, 
yn-2'®h-2^ ■ ^11-3’ %.- 3 ^] 


Corrector ; 

,(« 




[9f(Xji+l’ ^n+1’ ^+1^ 


“ 5 ffXn _ i > yn_i> ■*■ ^^^n-2’ ^11-2’ ^-2^^ 
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APPENDIX n 

CMDU (RAPP^ REACIQR DATA 

The iollowing fiajasthan Atomic Power Plant data of a CAKDU 
tj^po reactor is used here for the required calculations. 

Each fuel bundle consists of 19 fuel rods of Natural 

UO,. 

(1) Reactor inlet temperature T^ = 480*^ 

(2) Reactor outlet temperature T 2 = 560^' 

(3) i-ican D 2 O temperature T = 515°P 

(4) ToBiporature at the centre of the core T^^ (calculated) 

= 580 °!' 

(5) B 2 O temperature at low power = 504^ 

(6) Specific heat at constant pressure of the DpO coolant G 

= 1.7744 Btu/lb-E° 

(7) Dynamic viscosity of the D 2 O coolant p = 0.2838 Ib/ft-hr 

(8) Thermal conductivity of D^O coolant K = 0,3066 - g 

hr-E -ft 

(9) Reactor chmmol length 2i = 16.404 ft 

(10) Moan DgO density f = 52,91 Ib/ft^ 

(11) Mean DgO density in the riser at outlet temperature 

(calculated) pj, = 52.91 

(12) Moan DgO density in the downcomer at the inlet 

temporature = 57. 29 Ibm/f t^ 

( 13 ) Total system losses (calculated) = 1,413 Ibf/in 

2 

(14) Driving pressure (calculated) = 0.1827 Ibf/in 
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(15) Pump halftime for a tj^pical recirculation pump, 

l/cu)o = 5 sec, 

(16) Height of tho riser (hotleg) and downcomer (coldleg) 

of the U-tubc heat exchanger = 6 ft 

(17) Prandtl number of pr = pC^/K = 1.089 

(18) Roinolds nimabor of Ho = fvD/p = 5.68 x 10^ 

(ior the coolant velocity 27.2 ft/sec) 

(19) Husselt number (for the coolant velocity 27.2 ft/sec) 

= Hu = 654 

2C) Heat transfer coefficient (calculated for the given 
Nusselt number) h = Hul./D = 9300 Btu/hr-ft^-P 

(21) Hydraulic diamotor D = 0.0216 ft 

(22) Steady state Grashoff number basing on the diameter 

gpd^ - 

ol tho fuel rod (calculated) Gr^ = ^2 

= 0.4262 X 10® 

(23) Kinematic viscosity “))= ^ = 5.3638 x 10 ft /hr 

r 

(24) Clad surface temperature = 573. 8°F 

-4 

(25) Thermal expansion coefficient p = 4.0 x 10 P 

(26) Volumetric thermal source strength at the centre 

of tho core = 1.648 Btu/ft^-sec 

(27) Diffusivity a ~ 0.0032 ft /sec. 






Fig. 5-1 


G^ncrcil flow chQrt for moin 
prograrti applied to all the 

eq'U.dtiohS;';: 





Fig. 5'2 Flow chart of Integration 
Subroutine Adams 
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